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ABSTRACT 



In order to do relativistic gravimetry one needs to define 
a system of null coordinates for a given constellation of 
satellites. We present here three methods in order to find 
the null coordinates of an event in a Schwarzschild ge- 
ometry. We implement these three methods for the weak 
gravitational field of the Earth, compare their precision 
and time of computation. 

Key words: General relativity; relativistic coordinates; 
relativistic gravimetry. 



1. INTRODUCTION 



The classical concept of positioning system for a Global 
navigation satellite system (GNSS) would work ideally 
if all satellites and the receiver were at rest in an iner- 
tial reference frame. But at the level of precision needed 
by a GNSS, one has to consider curvature and relativis- 
tic inertial effects o f spacetime, which are far from being 
negligible llAsh03ll . There are at least two very different 
ways of dealing with relativistic effects. 

One way is to try to preserve the Newtonian conception 
of absolute time and space, by adding corrections com- 
ing from general relativity. This leads to numerous cor- 
rections that depend on the kinematics of the positioning 
system (see llAsh03fl for a detailed description). How- 
ever, the natural evolution of GNSS is to become more 
and more accurate with the help of very stable clocks, 
and autonomous (no need for ground stations) with cross 
links between the satellites. Galileo will embark hydro- 
gen maser clocks with a drift of about 1 ns after one day. 
But state of the art atomic clocks are far more stable than 
this (e.g. a drift of 26 ps/day for the Cesium clock Pharao, 
and only 0.3 ps/day for optical clocks). At this low level 
of uncert ainty a lot of supplementary corrections have to 
be added llLT02ll . 

Another way to define a positioning system is to abandon 
the Newtonian concept of absolute space and time, which 
is known to be a classical approximation, and to define a 
relativistic positioning system , with the so-called emis - 
sion, GPS or null coordinates IIBGHO02L ICPOfil lRov02H . 
A project called SYPOR, "SYstme de POsitionnement 



Relati viste", has bee n proposed by B. Coll and collabo- 
rators ||Co103L |Pas07fl . The relativistic positioning system 
is composed of four satellites broadcasting their proper 
times by means of electromagnetic signals. The coor- 
dinates of the receiver are simply the four proper times 
{n, T2, Ts, T4} received from the four satellites. Then the 
receiver knows his trajectory in these null coordinates. 
Moreover, if each satellite broadcasts its proper time to 
the other satellites, the receiver knows also the trajecto- 
ries of the satellites, and the system is auto-locating. 

The null coordinates have very attractive properties. First 
of all, they are covariant quantities, although dependent 
on the set of satellites one chooses and on their trajec- 
tories llLac06h . The set of satellites constitutes a pri- 
mary reference system, with no need to define a terrestrial 
reference frame; so there is no need to track the satel- 
lites with ground stations or to synchronize the clocks. 
There is no need also for relativistic corrections, as rel- 
ativity is already included in the definition of the posi- 
tioning system. If needed, they can be related to more 
usua l terrestria l coordina te system. Coll and collabora- 
tors IICFM06bL ICFM06aH studied such relativistic posi- 
tioning systems in the case of a two-dimensional space- 
time, for geodesic emitters in a Minkowski plane and for 
static emitters in the Schwarzschild plane. The relativistic 
positioning system has been studied in the vicinity of the 
earth, performing calculati ons at f i rst ord er curvature in a 
Schwarzschild spacetime [BahOl, RT08]. On one hand, 
such an approach at first order does not take advantage of 
the full meaning of a relativistic positioning system, and 
suppose that an underlying (non physically sounded) co- 
ordinate system is predefined; on the other hand, it has the 
advantage to give a model of the spacetime geometry in 
the vicinity of the earth, to which the data of a relativistic 
positioning system can be compared. 

The next generation of GNSS will have cross-link capa- 
bilities. Each satellite will broadcast the proper time of 
the other satellites in view, as well as their proper time. 
With these informations, one could in principle map the 
spacetime geometry in the vicinity of the constella tion of 
satellites by solving an inverse problem flTKPC09H . This 
proposal is at the origin of this study. Here, the use of rel- 
ativistic coordinates is not proposed to enhance the pre- 
cision of a positioning system, but to achieve a scientific 
objective, which is the mapping of the gravitional field in 
the vicinity of a GNSS. In the inverse problem, the basic 
unknowns are the components of the spacetime metric 



in the null coordinate system. A primary constellation, 
consisting of four satellites, defines the null coordinate 
system; the secondary constellation is constituted by the 
satellites that do not contribute to define the null coor- 
dinates. The data of the problem are mainly the proper 
times of the satellites broadcasted to the other satellites. 
These signals are electromagnetic waves, so their trajec- 
tories are null geodesies. Knowing the trajectory of the 
satellites and the proper time of emission of a signal, one 
can solve the geodesic equation and compute a proper 
time of arrival; by comparing the computed proper time 
with the observed one we can optimize our knowledge 
of the spacetime metric. In case the satellites are not in 
free-fall, additional data coming from an accelerometer, 
a gyrometer or a gradiometer could be used. Then the re- 
constituted geometry can be compared to a parametrized 
relativistic model of the spacetime geometry around the 
earth. 

As the existing GNSS are not autonomous, the required 
data to do a mapping of the gravitational field are missing. 
We thus propose as a project going beyond this article 
to simulate these data, and to apply an inverse problem 
in order to recover the initial metric. In the first part of 
this article we introduce basic results of null coordinates 
in a flat spacetime. In the second part we present three 
methods in order to solve numerically the time transfer 
problem in a Schwarzschild spacet ime: the first method 
relies on the world function and requires the im- 

plementation of a numerical integration algorithm and a 
shooting method; the second method relies on an ana- 
lytic solution of the geodesic equations in terms of elliptic 
functions [CK05]; and the third method relies on series 
in power of G - t he gra vitational constant - of the time 
transfer function llTL08ll . In the third part of this paper 
we apply these three methods to find the null coordinates 
of a stationary point in a Schwarzschild spacetime, we 
compare them in terms of precision and time of compu- 
tation. 



2. BASIC RESULTS IN FLAT SPACETIME 



The null coordinates in a four dimensional flat space- 
time has been studied by several authors ( 1BGRT081 
ICFM091 iBGHOOl iRovOa lRT08ln . We summarize here 
the main results. Let call x = (x 1 , x 2 , x 3 , x A ) the usual 
minkowskian coordinates, where x 4 = ct. The space- 
time metric in a flat spacetime is the Minkowski metric 
r] a (3 = diag[-l, -1, -1, 1] and 



ds 2 = r] a pdx a dx^ 



(1) 



We consider four satellites A = {1..4} that constitute 
a relativistic positioning system. We assume that each 
satellite has on board a perfect clock that delivers its 
proper time t a . Then the geodesic equation for the tra- 
jectory xa(t a ) of satellite A, parametrized by its proper 
time, is 



d 2 x A 

d(T A ) 2 



= 0. 



(2) 



The solution of which is: 

x a {t a ) = U A r A + S° A 



A • 



(3) 



where U A = dx A /dr A is a normal constant vectofl To 
simplify further we assume that S A = 0. 



Satellite 1 



Satellite 2 




X 



Figure 1. Illustration of the problem in two dimensions: in order to 
find the null coordinates of point P, one has to find the intersection of 
its past cone with the satellite trajectories. 



Let P, the point of coordinates xp, be the receiver lo- 
cation. To find its null coordinates we need to find the 
null geodesies linking this point to the four satellite tra- 
jectories (see Fig ure [p. To do this one can use the world 
function !lSvn60tl (see appendicelAlfor a definition). The 
condition for the signal emitted by the four satellites to 
meet the receiver location is 

n(x A (r A ), x P )=0, x\(t a ) < x P . (4) 

In flat spacetime the world function is 



n(x P , x A ) = -p a p (x P - x A ) (x P - x A ^j . (5) 

Then, we find 

\x P \ 2 -2t a (x p ,U a ) + (t a ) 2 = 0, (6) 

where we define (a,b) = i] a f3a a bP and |a| 2 = (a, a). 
We choose the solution of this equation which satisfies 
the inequality in Q. It defines the null coordinates: 



x P ,U A ) - \J (x p ,Ua) - \xp\ 



(7) 



This relation is true for any point P of spacetime: then 
the coordinate transformations r = t(x) (where r = 
(t 1 , t 2 , t 3 , t 4 )) are known everywhere, and the compo- 
nents of the contravariant metric in the null coordinates 
are 



Q T A Q T B 



dx a dx? 



(8) 



'Normal means that r\ a pU c pJJ® A = 1. 



In order to find the components of the covariant metric 
gAB one can find the inverse metric, or solve the sys- 
tem ( H} with xp as the unknown. As Bini et al. empha- 
sized llBGRT08ll . it is equivalent to solve the system: 

(r 1 ) 2 + |:zp| 2 - 2T 1 (x P , U A ) = (9) 
(r 1 ) 2 - (r 1 ) 2 - 2 (r* - r 1 ) (x P , U A ) = (10) 

where i = 2. A. In the system of equations (l9l-(fTQt there 
is only one equation of second order, instead of four for 
the system ©. Then 

dx a dx 



9AB 



8t a 0t b ' 



(ID 



Alter nativel y, one can use the time transfer func- 
tion llTL08ll to find the null coordinates of point P (see 
appendice [A] for a definition). In a flat spacetime, the 
time delay function is trivial (eq. (l38l l) so that 



T 4 



■X i A(T A )=R AP (T A ), 



(12) 



which is equivalent to the equations ©. 



SOLVING THE TIME TRANSFER IN A 
SCHWARZSCHILD GEOMETRY 



3.2. Method 2 



The second method is to solve the time transfer equation 
using the analytic solution of the geodesic equations in 
terms of elliptic functions llCK05ll . We summarize here 
the main results. The differential equation for orbits of 
lightlike geodesies is 



du 
dA 



± v /a 2 -u 2 (l-u), 



(18) 



where u = rs/r, A is the true anomaly of the orbit and 
a is a constant of motion. The solution of this equation 
for type A orbits (ie. scattering orbit with both end points 
at infinity) can be expressed in terms of elliptic functions 
with the transformation u = t 2 (u2 — U3) + 113: 



u 2 — («2 — "3)cn 2 



F{xA\m) + 



A - A/ 



(19) 



where m = (1x2 — Ui)/{u\ — U3), n = 2/ ^/u\ — 113, 
and u\, 112 and 113 are solution of the cubic equation a 2 — 
x 2 (l — x) = 0; A a is the true anomaly of a point A lying 
on the trajectory, and \A is such that 



The spacetime metric written in the usual Schwarzschild 
coordinates x = (r, 9, </>, ct) is 

ds 2 = g a pdx a dx 



= A(r)dt 2 - A' L {r)dr 2 - r 2 dVL 2 



where 



A(r) 



1 



rs 



(13) 
(14) 



(15) 



dO 2 = (d6 2 + sm 2 (8)d(j) 2 ) and r s is the Schwarzschild 
radius. 

Let call xp the coordinates of the point of reception 
P, and Sa = xa(t a ) the worldline of satellite A, 
parametrized by its proper time t a . 



3.1. Method 1 

The first method to find the null coordinates r of point P 
uses the world function: one has to solve the system of 
equations @. This system here is equivalent to 



dA 



= 0, 



(16) 



where A is an affine parameter, chosen such that xl(0) = 

xa(t a ) and xl(1) = xp, () = d/dA, and xl(\) is 
solution of the geodesic equations: 



dA 



0, 



(17) 



where V is the covariant derivative and ul = dxi, /dA. 
This forms a two point boundary value pr oblem t hat can 
be solved using a shooting algorithm (see llSan07ll ). 



. 2 U A 

sin xa 



U3 



U2 - U3 



(20) 



where ua = rs/r a, with the radial coordinate of 
point A; cn is the usual elliptic function and F is the in- 
complete elliptic integral of the first kind. 

If the spatial coordinates of two different points A and 
P are given, then the parameter a of the null geodesic 
crossing these two points can be found by solving a non- 
linear equation on a, which can be found by using the 
Jacobi elliptic functions addition theorem and the general 
expression (O of the orbit (eq.(19) of llCK05ln . Once 
the parameter a of the null geodesic between the points 
A and P is known, one can calculate the time of flight 
Tf(xp,XA]a) of the photon between A and P by using 
the formula (25) of flCK05ll . If the arrival coordinate time 
xp is given, and that the emission point A lies on the 
given orbit xa(t a ) of satellite A, then one can find the 
null coordinates r of P by solving numerically 



x A (T A ) = T f (x P ,x A (T A );a). 



(21) 



3.3. Method 3 



The third and last method uses a post -Mink owskian ex- 
pansion of the time transfer function lTL08ll for a static 
and spherical body. The second order post-Minkowskian 
expansion of the time delay function between two points 



A and P (see appendice[A} is then (for general relativity): 



A(x' Al x' P ) = r s ln 



R' 



AP 



R' 



AP 



„2 R'aP 



15 arccos(n^ • n' P ) 



l + (n' A - n' P ) 



(22) 
(23) 



where x' are the quasi-Cartesian isotropic coordinates, 



(with the euclidean norm), r f 



r 'ap = W&a - f pll> n A = x'a/t'a' n'p = x'pl r 'p and 
r$ is the Schwarzschild radius. This function is indepen- 

dant of the time of reception and the time of emission, as 

it should be for a stationary spacetime. 

The transformation from the Schwarzschild coordinates 
to the quasi-Cartesian isotropic coordinates can be ob- 
tained with: t' = t, 8' = 9, <fr' = <fi and (up to the second 
order in (rg/r)) 



1 _ rs_ _ /rsV 
2r V4r/ 



(24) 



Then, for a given point P of coordinates x P , one obtains 
its null coordinates r by resolving the four equations: 



x P - x\t a ) = R ap + A (x(t a ), xp) 
where x{t a ) is the given orbit of satellite A. 



(25) 



Note that this method is only valid to second order in 
(rs/r). It is possible to obtain the higher order terms in 
equation ( 1221 but at the price of complicated calculations. 
On the contrary, the precision of the two first method is 
only limited by the time of computation of the numeri- 
cal method used to solve the problem. In the next part 
we compare the numerical implementation of these three 
methods to find the null coordinates of a stationary point 
in a Schwarzschild spacetime. 



THE NULL COORDINATES OF A STATION- 
ARY POINT IN SCHWARZSCHILD SPACE- 
TIME 



4.1. The constellation of satellites 



We want to study a constellation of four satellites in the 
plane 8 = ir/2, with circular orbits, and that define a 
null coordina te syste m. The geodesic equations for the 
satellites are llPap74ll : 



'0 




2r 3 



(26) 
(27) 
(28) 
(29) 



where () = d/dr, with t an affine parameter. Let t be 
the proper time, then we have the supplementary equation 
u a u a = 1, where u a = dx a /dr. The solution is: 



r = r 
t(r) = 



1 



3rs 
2r 



-.T + t 



(?r s 



2rl 



3rs_\ 
2r ) 



(30) 
(31) 

(32) 



where to and 0o ar e two constant of integration. No 
circular orbits ex ist for ro < 3r$/2. It is also well- 
known IMTW73H that for 3r s /2 < r < 3r 5 , circular 
orbits are unstable, and stable for tq > 3rs- One can 
notice that 



(33) 




2rl 



In order to compare the different methods sketched in the 
previous section, we need to solve the problem for only 
one satellite A of the constellation. Its initial conditions 
are ro, to an d </>o- To simplify, we choose to = and 
4>q = 0, and we note ro and t a the radius and the pa- 
rameter of the satellite trajectory. Let P be a point of 
given coordinates x P , and xa(t a ) the trajectory of the 
satellite. 



4.2. Numerical results 

The first algorithm uses a shooting method to solve the 
Two Point Boundary Value Problem (TPBVP). The TP- 
BVP is formulated as an initial value problem (IVP), 
where the initial point has to be guessed for the termi- 
nal point to satisfy the terminal conditions. The IVP can 
then be solved using a Newton-Raphson algorithm. To 
compute the Jacobian we use finite differences. An al- 
ternative can be computing the transition matrix, which 
would result in the integration of a 9 2 differential system. 
When using finite differences, the time step has to be se- 
lected carefully. One is always tempted to take the small- 
est step possible, although that does not imply a better 
approximation of the derivatives. One elegant approach 
is to have an accurate numerical value of the derivatives 
using the theory on holomorphic functions. Although, 
the functions have to be extended to the complex space, 
which is not always readily feasible. Another option is 
to use automatic differentiation to help to reduce the high 
burden of calculating the analytical derivatives. 

Methods 2 solves two algebraic equations involving Ja- 
cobi elliptic functions and elliptic integrals. The equa- 
tions are unfortunately not well suited numerically to be 
solved readily with N ewton-type solvers. We use the 
Brent method llBre73ll . which is a secant method algo- 
rithm, to solve the problem. Another option would be 
to expand the integrand to calculate the time of flight Tf 



(formula (23) of llCK05lO into a convergent series of ana- 
lytically integrable functions, which should be more suit- 
able and faster for weak gravitational fields llCK05ll . 

Method 3 does not require any integration either. We sim- 
ply need to solve one algebraic equation that happens to 
be a polynomial equation of second order. This equation 
can be solved efficiently with a Newton-type algorithm, 
where the Jacobian can be computed accurately by finite 
differences. 



In this study we used only double-precision calculation. 
To have the most accurate results and have a fair compar- 
ison between the methods, we should use multi-precision 
calculation for the three methods, but its implementation 
is not yet done for the first two methods. In a multi- 
precision implementation with a 128 bits description of 
the real number, we should expect at least 30 significant 
digits. 

We computed with the 3 methods the null coordinate t a 
of the point P of spatial Schwarzschild coordinates 

r P = 50000 km 
Bp = it/2 

4>p = o, 

for different arrival times tp. The initial conditions of 
satellite A are 

r Q = 42000 km 
to = 
0o = 0; 

its trajectory is given by equations (I30]l-(l32l and is sam- 
pled with about 100 points. We choose the Schwarzschild 
radius of the Earth. 

We use as initial guess, for all three methods, the solu- 
tion to the flat space problem. This reduces the compu- 
tation time and ensure a good convergence. The results 
are shown on figure [2] As the 3 methods provide similar 
results to the first digits, the graphs are superimposed and 
no difference is visible. 

To compare the methods, we take 4 points sampling the 
satellite orbits, compute the discrepancies between the 
solutions for each method, and evaluate the computation 
time. We take as reference value, the value returned by 
the method 1. The results are displayed on tableQ] 

The computational times (TC) figured on table Q] as ra- 
tio with respect to the computational time of method 1. 
These durations are obtained for double-precision num- 
bers, so not for the most accurate comparison we could 
do between the different methods. Assuming that an 
efficient multi-precision implementation of the methods 
would lead to similar results, we can infer that method 3 
is the fastest. However this method will give results pre- 
cise up to the second order in (rs/r), and we are inter- 
ested in a method that gives a result to any desired accu- 
racy. Methods 1 and 2 are equivalent regarding the preci- 
sion, but method 1 is hardly faster than method 2. This is 



Pt 


tp(s) 


t a (ref) 
Method 1 


Method 1/2 


Method 1/3 


1 


10° 


.9733148699 


9.8861e~ i;a 


7.801e- i;> 


2 


10 1 
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1.0181e" i3 
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10^ 


99.9732913262 


2.3521e" Y 


8.9951e-^ 


4 


10 a 


999.9710561425 


2.2591e- b 


7.12916- 11 


TC 




1 


1.25 


0.5 



Table 1. Numerical comparisons between the methods 



mainly because using a Brent algorithm is far from being 
efficient regarding computation time. 




Figure 2. Coordinate t a of the stationary point P as a 
function of the arrival time tp. 



5. CONCLUSION 



In this paper we recalled the basic properties of null co- 
ordinates in flat spacetime. We presented three meth- 
ods in order to find the null coordinates of an event in 
a Schwarzschild spacetime. We implemented these three 
methods and compared them in terms of time of computa- 
tion and of precision. The post-Minkowskian expansion 
(method 3) is faster than a direct integration (method 1) 
or than an evaluation of the elliptic functions (method 2). 
However using only double precision in the calculation 
leads to time of calculation that are of the same order. 
We still need to implement the three methods in multi- 
precision to compare them to the maximum accuracy of 
the third method. To go beyond this preliminary study, 
we want to use these results in order to simulate the data 
of a constellation of satellites in the context of relativistic 
gravimetry. 



A. DEFINITIONS 

The world function llSvn60ll Let P and Q two points 
of coordinates x p and xq, joined by a geodesic T : x = 



£(A) where A is an affine parameter; then the world func- 
tion O is 

1 r 1 

Q(x P ,x Q ) = - j g a0 {x)\ v U a UP&\, (34) 
where U = d£/dA and A is such that £(0) = xp and 

€(1) = *Q. 

The time transfer function dTL08ll We assume that 
spacetime is globally regular and without horizon, and 
that the coordinate system is chosen such that 344 > 
everywhere. The past null cone at a given point P of co- 
ordinates xp intersects the worldline Sa at one and only 
one point xa = (xA,ctA) (where xa = (x\,x A ,x A )). 
The difference tp — tA is the (coordinate) travel time of a 
light ray connecting the emission point xa and the recep- 
tion point xp. This quantity may be considered either as 
a function of the instant of reception t p and of xa, xp, or 
as a function of the instant of emission tA and of xa and 
xp. So it is possible to define two time transfer functions, 
%■ and T e with: 

tp -t A = T r (xA,tp,x P ) = T e (t A ,XA,x P ). (35) 

% is the reception time transfer function and T e the emis- 
sion time transfer function. These functions are distinct 
except in a stationary spacetime in which the coordinate 
system is chosen so that the metric does not depend on 
x A . In this case the subscript r and e can be omitted. 

We suppose now that the metric takes the form 

9ap = Va/3 + h a /3. (36) 

Then the reception time transfer function may be written 

as 

T r (xA,tp,x P ) = -Rap + -A r (x A ,tp,x P ), (37) 
c c 

where Rap = \\xa — xp\\ (with the euclidean norm). 
The function A r /c is called the reception time delay 
function. 

In fiat spacetime it is obvious that 

A(x A ,xp) = i). (38) 
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